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In a recent work, M. Troyer, F. Alet and S. Wessel Q] proposed a way to extend histogram methods 
to quantum systems in the World Line Quantum Monte Carlo (WLQMC) formulation. The strategy, 
also proposed in 0, allows to compute quantum averages on a narrow temperature range from a 
single Monte Carlo run at fixed temperature. This is achieved by fixing N, the number of temporal 
divisions in the Trotter-Suzuki expansion of WLQMC, and by changing £=l/(JVfcjT). In this work 
we apply this strategy to construct a single histogram Monte Carlo method for a canonical ensemble 
of one- dimensional quantum harmonic oscillators and we explore its accuracy limits. We obtain 
that fixing N imposses a limit of minimal temperature to the properly performance of the method, 
which is T m i, I = 1.9(2)7V~ 0,80l ' 6 ' ) in our example. This limit is a consequence of the fact that the 
Trotter-Suzuki expansion fails for large e values, and, therefore, should be taken into account in all 
applications of this histogram method for quantum systems. 
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INTRODUCTION 

Histogram Monte Carlo methods 0- are pow- 
erful and robust techniques to investigate classical sys- 
tems, both with discrete or continuous degrees of free- 
dom. Such methods have many advantages over tradi- 
tional Monte Carlo ones. For example, one can obtain 
mean values of thermodynamical quantities over a broad 
range of temperatures by using data from a single run. 
Also it is possible to compute thermodynamic potentials, 
like the free energy, that are not accessible by other me- 
thods. All these methods bases on defining a density of 
states, let say g(E), i.e. the number or states with the 
same energy. Once obtained, g(E) allows to compute av- 
erages like the mean energy and the specific heat at any 
desired temperature. Therefore, the problem of averaging 
over all states is transformed into a problem of counting 
to find g(E). Different histogram methods differ in the 
way they obtain the density of states. 

There are histograms methods for analysis and his- 
tograms methods for sampling. Some examples of the 
first ones are the Single Histogram Method (SHM), first 
proposed by Z. W. Salsburg in 1959 0] and populari- 
zed by A. M. Ferrenberg and R. H. Swendsen in 1988 @; 
the Broad Histogram Method (BHM) formulated in 1996 
by P. C. M. de Olivcira, T. J. P. Pcnna and H. J. Her- 
mann 0, and the Transition Matrix Monte Carlo me- 
thod (TMMC) proposed by J. S. Wang, T. K. Tay and 
R. H. Swendsen in 1998 [1 Q, which reduces to BHM 
when a single energy jump is chosen. In the second class 
we find, for example, the Multicanonical Method, pro- 
posed by B. A. Berg^ (also formulated by J. Lee as 
the Entropic Sampling metho dllll L the Flat-Histogram 
method of Wang and Lee 0, Eland the Wang-Landau 
method [3| . These methods take samples uniformly dis- 



tributed on the energy axis. 

Recently, M. Troyer, F. Alet and S. Wessel proposed 
a way to extend histogram Monte Carlo methods to ca- 
nonical emsembles of quantum systems, in the world line 
quantum Monte Carlo (WLQMC) formulation In 
contrast withprevious proponsals 0, ^| , this idea, also 
proposed in [2j resembles the full power of the clasical his- 
togram methods, since it is able, for instance, to calculate 
mean values on a narrow temperature range from a sin- 
gle Monte Carlo run at a fixed temperature (SHM). The 
idea bases on the observation that interactions in time 
and space in the equivalent classical system of WLQMC 
are weighted by the temperature in different ways, as we 
will illustrate below. Thus, it is useful to define two new 
quantities, k\ and k%, which summarize spatial and tem- 
poral interactions, and also a density of states g{ki,k,2) 
on them. Such density can be computed by using any 
histogram method and allows to obtain averages at dif- 
ferent temperatures from a single run. This strategy has 
been used by Troyer, Alet and Wessel to investigate the 
three-dimensional quantum Hciscnbcrg model [l|, with 
excellent results. 

In this work we illustrate this strategy by studying a 
much more simple system: a canonical ensemble of non- 
coupled onc-dimcnsional quantum harmonic oscillators. 
By doing so, interesting aspects of the method like its 
range of accuracy and its dependence on the simulation 
parameters will be explored. The article is organized as 
follows. First, the World Line Quantum Monte Carlo me- 
thod is revised. Next, the quantities k\, ki and g(k\,k2) 
for a quantum single-particle system are defined, and it 
is shown how this density of states allows to compute 
expectation values for the kinetic, potential and total en- 
ergy operators at any desired temperature. Then, we im- 
plement the SHM for the one-dimensional quantum har- 
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monic oscillator, and we investigate how the temperature 
range of accuracy is affected by the number of divisions 
on the axis of imaginary time. Finally we summarize our 
main results and conclusions. 



THE WORLD LINE QUANTUM MONTE CARLO 
METHOD 



The World Line Quantum Monte Carlo method 
(WLQMC) was formulated in 1982 by J. E. Hirsch, 
R. L. Sugar, D. J. Scalapino and R. Blankenbecler [lgj . 
It has been used to study one particle systems as well as 
many particle systems, like interacting bosons and sys- 
tem of strong correlated fermions. The WLQMC em- 
ploys path integrals to change the partition function of 
a quantum system to the partition function of an equi- 
valent classical system in one higher dimension 0|. To 
illustrate the method, let us consider a quantum system 
with a separable Hamiltonian 



H = T(P) + V(X) = — 



V(X), 



(1) 



where V(X) and T = P 2 /2m are the potential and ki- 
netic energy operators respectively. In the canonical en- 
semble, the partition function is 



Z = Tr[e-P H ] 



(2) 



where — is the imaginary time, H is given by 
and is the Boltzmann constant. 

Since X and P do not commute, the exponential in 
@ cannot be separated into a product of exponentials. 
Therefore, the imaginary time is written as = Ne 7 
and a Trotter-Suzuki approximation is used |2ll |22J, l23j | 
to write the partition function J2J) as 

Z= lim Ztr= lim Tr[( e - ef e -6 ^* 5 )*]. (3) 

N—>oc N—^oo 

To compute lf!T|. the trace is expressed as a sum over a 
complete set of position eigenstates, and a new complete- 
ness relation is introduced between each of the N factors 
in that expression, obtaining 



= [^-J / Vx exp[-£ff(s)], 



(4) 



JV-l 



i N ~ 1 / _ \ 2 
H{ x ) = Y,V{x l ) + \mY J [^ J ^ 1 -) , (5) 

where T>x = dx^dxi ■ ■ ■ cLxn-i- The partition function 
for the quantum system is now expressed as a multidi- 
mensional integral over classical variables. The expres- 
sion J^J) can be interpreted as the energy of a chain of 
N classical oscillators with potential energy V{x) cou- 
pled through springs of constant k = m/e 2 . In this case, 



x = {xo,Xi, . . . ,xn-i} will denote a specific configura- 
tion of this classical system. 

By following the same procedure described above, one 
obtains 



» JVx V(x )exp[-eH(x) 

[ Hr JVx exp[-eH{x)} 



(G) 



(f)tr = 



2e 



m J Vx (xp - X!) 2 cxp[-eH(x)} 
2e 2 JVx exp[-eH(x)] ' { ' 

for the expectation values of the potential and kinetic 
energy operators respectively, where xq and x\ can be 
replaced by any pair (xi,Xi+i) of neighboring classical 
oscillators. For the total energy we have (H)tr — (T)tr + 
(V(X)) tr . 

Expressions © and Q show that the expectation va- 
lues for quantum operators are now expressed as multi- 
dimensional integrals over the configuration space of the 
classical system. This integral can be computed by us- 
ing standard Monte Carlo methods. For example, we 
can use the Metropolis algorithm to compute the to- 
tal energy average (H)t r as follows. First, an oscilla- 
tor i in configuration x is chosen at random and its 
position Xi is changed to a new value x\ + r, where 
r is a uniformly distributed random number in [—1,1]. 
Next, this new configuration x 1 is accepted with rate 
A(x' | x) = min{l, e -e '- ff ' a! ^ H< - x ^}. Finally, energy av- 
erages are computed as arithmetic means over the sample 
set. 

The temperature can be fixed by varying either N or 
e, since = 1/ksT = Ne. When N is fixed and e varies 
to fix the temperature, errors in the Trotter-Suzuki ap- 
proximation are greater at low temperatures. 



HISTOGRAM METHODS FOR QUANTUM 
SYSTEMS IN THE WLQMC FORMULATION 

For a classical system in the canonical ensemble, the 
average energy at temperature T is given by 



(Eh 



E E 9(E)Ee 
E E 9(E)e- 



(8) 



where g(E) is the number of configurations with the same 
energy E. If g(E) is known, one can use|Hlto obtain (E)t 
as a continuous function of temperature. So, the problem 
of summing over all states has been transformed into a 
problem of counting to find g(E). 

In the classical SHM, samples are taken with proba- 
bility proportional to e -Bx / fcsT (canonical distribution) 
at fixed temperature T, and a histogram of visits V(E) 
is accumulated as the number of samples with energy 
E. From such histogram, the density of states g(E) is 
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approximated as 



1(E) 



exp[E/k B T], 



(9) 



where J\f is the total number of samples. By replacing 
in JSJ), we can obtain averages at temperatures T' ^ T. 

To extend histogram methods to quantum systems, one 
has to define a density of states. As proposed in [llQ, 
this is achieved by defining two quantities in the WLQMC 
classical system, which resumes the interactions in space 
and time directions: 



N-l 



N-l 



h = J2 v ( x ^ k ^J2 & - x '+i) 2 ' ( 10 ) 



1=0 



1=0 



with periodic boundary conditions. The partition func- 
tion and the averages © and can be written in 
terms of k\ and k 2 as 

/ m \ N / 2 r 

Z tr = (^) / dk lC lk 2 g{k u k 2 )e-* H ^< k *\ (11) 



(V)t 



1 \dk Y dk 2 g(ki,k 2 )k 1 e ~ eH ^ k ^ 
N Jdkidk 2 g(ki,k 2 )e- eH ( k ^) ' 



(12) 



m J dk x dk 2 g{k u k 2 )k 2 e-' LH ^ klM ^ 



2e 2Ne 2 J dk x dk 2 g(k 1 ,k 2 )e- eH ^ k ^ ' 

^ (13) 

where H(ki,k 2 ) = k\ + mk 2 /2e 2 , and g{k\,k 2 ) is the 
number of configurations x = {xq, x±, . . . , xn—i} of the 
classical system with same values of k\ and k 2 for fixed 
N. 

Expressions (|1 ip — (|13f> are continuous equivalent to 
those of classical systems (JBJ • We can use any histogram 
method to compute g(k\,k 2 ). Let us use the SHM to 
illustrate this point. First, we fix the simulation tem- 
perature T s i m = l/Ne S i m by choosing N and the initial 
e s im- Then, samples are taken with probabilities pro- 
portional to e~ eH ^ by using the Metropolis algorithm 
described above, and an histogram of visits V(k\,k 2 ) is 
cumulated as the number of samples with values k\ and 
k 2 . Next, the density of states g(k\ 1 k 2 ) is estimated as 

g(fc 1; fc 2 )- V( y 2) exp[ £ g(fc 1 ,fc 2 )]. (14) 

Once g(k\,k 2 ) is computed, averages at different tempe- 
ratures are obtained by changing e in equations l(T^|) - ((T3^l . 
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Figure 1: Average potential energy for the one-dimensional 
quantum harmonic oscillator at various simulation tempera- 
tures. 
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Figure 2: Average kinetic energy for the one-dimensional 
quantum harmonic oscillator at various simulation tempera- 
tures. 



For this system, k\ takes the form k\ = ^mu> 2 J^ILo* x f> 
with m the mass and u> the frequenc y. T he exact expres- 
sions for the expectation values are |20|, |24| 



(H) = 2(T) = 2(V) = fiu(- 



1 



(15) 



Figures [5] and [3J shows the expectation values of 
the potential, kinetic and total energy, at three different 
simulation temperatures, T s , m . The parameters for each 
simulation are shown in table [I] In all cases, the total 
number of samples was 60000, ks = 1 , 30000 mess were 
discarded before sampling, and 10 mess were performed 
between successive samples. To make statistic, the whole 
procedure was run ten times with different seeds of the 
random number generator. 

Errors in the Trotter-Suzuki aproximation are propor- 
tional to e 2 01 ■ Since the SHM fixes N and changes the 



THE SHM FOR THE ONE-DIMENSIONAL 
QUANTUM HARMONIC OSCILLATOR 

The system we chose is a canonical ensemble of non- 
coupled one-dimensional quantum harmonic oscillators. 



Table I: Simulation Parameters 

N es^rn Tsirn C PUtimejs) 

40 0.25 0.10 412.294 

30 0.09 0.37 312.607 

10 0.10 1.00 117.532 
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Figure 3: Average total energy for the one-dimensional quan- 
tum harmonic oscillator at various simulation temperatures. 
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Figure 4: T min as a function of N for both SHM and WLQMC 
in log — log scale. 



temperature by variying e, the errors will be greater for 
lower temperatures. There is a minimum temperature, 
Tmin, where the simulation results includes the theorical 
ones between the error bars. Such a T m i n is a function 
of N. By looking at the total energy, we found for the 
SHM 

T mm = 1.9(2)iV- - 80 ( 6 > . (16) 
Following the same procedure for the WLQMC we found 

T mm = 2.0(6)iV- a81(6) . (17) 

This two results are very similar and are plotted in 
Figure They actually show that the minimum tem- 
perature in SHM is a consequence of the errors in the 
Trotter-Suzuki aproximation. This limit is important 
for the correct performance of the method and should be 
taken into account for all applications of histogram me- 
thods for quantum systems in the WLQMC formulation. 

Increasing N or decreasing e gives exacter results at low 
temperatures, but at price of longer CPU times, just be- 
cause correlation times increases with decreasing e [lflj . 
Nevertheless, the use of SHM spares large amounts of 
CPU times. For example, it takes 118 sec on a Pentium 
4 at 1.6 GHz to compute (H)t over the whole tempe- 
rature range 0.7 < T < 1.7, at simulation temperature 
T s im = 1.0, in contrast with the 35 sec of a single point 
by WLQMC. 

CONCLUSIONS 

In this work we have employed the extension of his- 
togram methods for quantum systems in the WLQMC 
formulation proposed in 0, 0] to the case of a canoni- 
cal ensemble of one-dimensional quantum harmonic os- 
cillators. By doing so, we have explored the ranges of 
accuracy of the method as a function of the number of 
divisions N in time direction. Actually, since the me- 
thod fixes N and varies e to change T, the corrections 



in the Trotter-Suzuki expansion fixes a minimum tem- 
perature for the correct performance of the SHM. We 
obtained that the method is accurate for temperatures 
above T mm =1.9(2)A-°- 80(6) . That is corroborated by 
direct WLQMC simulations, cciuationll7l This T m i n can 
be lowered by increasing N, but at cost of larger CPU 
times, because correlation times grow with decreasing e. 
Perhaps SHM needs smaller CPU times than WLQMC. 

The definitions of k\ and ki and the density of states 
g(ki,k2), expressions (fTU)l and (fl"4"jl . allow to apply any 
histogram Monte Carlo method to quantum systems in 
the WLQMC formulation, with all the power and advan- 
tages that histogram methods have deserve for classical 
ones. However, the restrictions founded here applies to 
all of them, because they relies on the fail of the Trotter- 
Suzuki expansion for lower e values, and should be always 
taken into account. Further analysis, like the theoretical 
study of the errors in this method, are interesting areas 
of future work. 

The strategy of extending histogram methods to quan- 
tum systems in the WLQMC formulation is completely 
general, and is able to give accurate results for many 
kinds of systems. We expect that our work will be useful 
to the right use of the method for many systems in the 
future. 
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